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I. INTRODUCTION 


Bi-orthogonality relationships for porous media have been developed by Pro- 
fessors Scandrett and Frenzen [Ref. 1]. They showed that Fraser’s [Ref. 2] bi- 
orthogonality relationship for the propagation of Rayleigh-Lamb modes in a plate 
with traction-free surfaces is a special case of the porous case in which the porous slab 
has zero porosity. Bi-orthogonality relationships can be used to solve scattering prob- 
lems from interfaces of discontinuity between vertically stratified porous/elastic/fluid 
media. In this thesis, we will examine the bi-orthogonality relationship for the elastic- 
fluid case. 

Before the bi-orthogonality relationship can be used, the eigenvalues and eigen- 
functions corresponding to the boundary value problem must be found. In this thesis, 
we will attempt to find these numerically. The question we ask ourselves is which 
numerical approach is the best. One method is to determine Rayleigh-Lamb modes 
[Ref. 3] in the elastic solid and try to perturb those solutions to obtain modes for 
the fluid loaded layer. The problem with this method is that modes may be missed. 
Another approach is to use the Rayleigh-Ritz method [Ref. 4], which is based upon 
the calculus of variations. It uses an expansion of elementary functions to “minimize” 
a functional to the corresponding boundary value problem. The functions must be 
carefully chosen due to the complexity of the resulting coefhcient matrix which will 
have terms involving the unknown eigenvalues. 

The method chosen in this thesis is the finite difference method. This method 
is very useful when numerically modeling boundary value problems involving deriva- 
tives. Using the finite difference approximations, we produce a discrete linear operator 
which approximates the differential operator of the boundary value problem. We can 
then use MATLAB [Ref. 5] to solve for eigenvalues of the linear operator which 
are theoretically close to the eigenvalues of the differential operator. We will also 


use additional numerical methods, to be discussed later, in an attempt to refine our 


approximate eigenvalues. 

Chapter II will provide a background in acoustical theory and will introduce 
terms that the reader will need to be familiar with for further chapters. Also, the 
concept of bi-orthogonality is explained in detail and a simple example is presented. 
Finally, a description of numerical techniques and methods that are used throughout 
the thesis will be given. 

In Chapter III, background on the elastic-fluid problem is introduced. Detailed 
analytical solutions to several elastic half space problems are provided to establish the 
basis of the numerical solution. The fluid loaded plate case will be the case studied 
in the numerical portion of this thesis. 

Chapter IV examines the numerical difference approximations used in the 
problem. The modal wave speeds are presented and compared to the “refined” wave 
speeds. Comparisons will also be made between the analytical and numerical eigen- 
functions/vectors of the fluid loaded elastic plate. 

Finally, Chapter V will provide conclusions and discuss areas of further re- 


search. 


ite THEORY AND METHODS 


A. ACOUSTIC THEORY 
1. Basic Acoustics and the Wave Equation 


Acoustics is generally associated with sound. However, the definition of acous- 
tics is the generation, transmission, and reception of energy that is in the form of 
vibrational waves. Audible waves are between 20 and 20,000 Hz where 1Hz = 
lcycle/second. The simple oscillator[Ref. 6], or basic spring mass problem, is the 
most common example used to demonstrate basic acoustics. The simple oscillator is 


governed by the equation 


he Te ane 
qe (11.1) 
with a general solution of 
xz = Acos(wot) + Bsin(wot) (11.2) 


where w is the angular frequency in radians per second (rads/sec). Since there are 
27 radians per cycle, we define frequency as 
Ww 
=—. IT.3 
j= (11.3) 
In addition, we note that the period T is equal to a If the initial velocity and 
displacement are known, we can obtain an exact solution to Equation (II.1). 

The fluid-elastic problem that is the central problem of this thesis is a two- 
dimensional problem, and thus the fluid and elastic potential functions are governed 
by the two-dimensional wave equation 

1 0? 


V*6 = ——_ I1.4 
Chae ee 


We assume an incident wave is travelling in the positive z-direction, which strikes an 
interface generating reflected and transmitted waves. A plane time harmonic wave 


propagating in the postive z-direction has the form 


® = plye rhs) (11.5) 


where ky = —. Substituting Equation(II.5) into Equation(II.4), one finds the ampli- 
C] 
tude function must satisfy the relation 


o'(y) + (ki — k*)y(y) = 0 (11.6) 


which is the canonical form of the boundary value problem to be solved. 


2.  Stress-Strain Relationships For Elastic Solids 
Elastic solids have the property that when external pressures are removed, the 
solid returns to its natural state. Strain is a measure of displacement while stress is 


a measure of force per area. We define the strain as 


1 
ij = 5 (ui + U3.) (11.7) 


where u is a displacement vector and tensor notation has been adopted. The term 
u;,; means the partial derivative with respect to z;, or for example, u),2 corresponds 


Ou ae 
to ——. The strain, €;;, corresponds to elements of a matrix. Note that ¢;; 1s dimen- 


Ox, 
sionless. Stress, 7;;, has the dimensions of force per area. For small strain, ¢, the 
stress, T, is proportional to strain. This relationship is known as Hooke’s Law [Ref. 


7], which can be written for an isotropic homogeneous solid as 


Tj = AERO; -- Yes ae (11.8) 


where A and pw, known as Lamé’s elastic constants, have the dimension of pressure, 


and 4;; is the well-known Kronecker delta defined as follows: 


i) 
Os = (11.9) 
0 273 
The displacement vectors can be represented in terms of potential functions. 


The displacement equation of motion for an isotropic homogeneous elastic solid [Ref. 


7] can be written as 


uV*ut (A+ p)VV-u= pi (1.10) 


We consider a decomposition of the displacement vector in the form 


u=Vy+Vx wo (11.11) 


It can be shown that the displacement vector u satisfies Equation (II.10) if 





Le. 
Vp = 9 
ob 
Le) 
ce 
Vp=-z¥ 
oR 
A+2 
where c? = a and G = F iRef. 7]. The question of the completeness of 
p p 


the solution in Equation (II.11) is guaranteed by the following theorem given in 
Achenbach’s book on elastic solids[Ref. 7]. 


Completeness Theorem: Let 


meee y x 1) 
fec(y xT) 


and satisfy the equation 
uy a + (A += 2)VV 2 + pf = pu 


in a region of space V and in a closed time interval 7’. Also let 


f =, VF ie ee Gr. 


Then there exists a scalar function y(z,t) and a vector-valued function o(Z,t) in 


terms of which u(Z,t) is represented by 
Ui VS Ne aa) 


with 


and where y(z,t) and w(z,t) satisfy the inhomogeneous wave equations 


L 
Wek Tee Tarra 2 
CL 


| = 
Vp+G= sy 
er 


The proof of this theorem can be found in [Ref. 7]. 

Representation of the displacement vector u in terms of potentials will become 
important in the remainder of this thesis. We restrict our analysis by assuming the 
strain in the z direction is zero, which corresponds to the two dimensional plane strain 


case. This allows us to specialize our decomposition into the following 


U= Yr t py 
V = Yy — vz 


(11.13) 


where the vector valued displacement potential is reduced to a single nonzero scalar 


component. Our stress equation follows from Equation(II.8) as 


Toy = p(y + ve) 
Ty = A(us + Vy) + 2pry 


(11.14) 


where, as before, the potentials y and w satisfy Equations (II.12) 


For an ideal fluid the stress in any direction equals the pressure stress of the 


fluid, or in equation form 


Trr = Tyy = —Pp (Je) 


B. BI-ORTHOGONALITY 

We say that two functions y(z) and w(x) are orthogonal over an interval 
Se [ y(z)v(x)dx = 0. This is analogous to the dot product of two 
orthogonal iis equaling zero. The orthogonality relationships[Ref. 8] for sines 


and cosines are as follows 














L 0 
I sin sin dz = oe (11.16) 
0 L L a ee 
0 ee ie 
L 
| cos—~cos——da =), 2 n=m#0 (II.17) 
aa. a) 
/ a ce (11.18) 
n— — 
: 7, cos > dx 


This orthogonality relationship allows us to determine values of coefficients when 
arbitrary functions are expanded in terms of these functions. For example, if f(z) 


can be written in the form 


TL 


AC » B, sin 7 (11.19) 


then we can apply the orthogonality of sines over the interval (0, L] to obtain 











L oe) ne 
i‘ f(z)sin—~de = a B, | sin sin da (11.20) 


which leads to 





L NTL 
f(x)sin dx 
z= AES (11.21) 


or the more familiar form 


NTL 


2 
B, == | f(2)sin de (11.22) 





The analysis is similar for cosines. Note that we have formally interchanged summa- 
tion and integration in the above steps which is valid provided the series and integrals 
converge. 

Bi-orthogonality relationships generalize the orthogonal relationships of clas- 
sical eigenfunction expansions. In [Ref. 1], these relationships were derived for waves 
propagating at a free surface and along interfaces between porous and elastic media. 
The general bi-orthogonality relationship for a two dimensional porous layer which is 


fluid loaded is: 


0 hy 
- eee a ee = Si.0) =| paenag “Ie | [(—p4U 7] pura y == (I) (Tie23} 
—fe2 


where the fluid layer and porous layer thicknesses are h; and hz respectively, and 
s = —(p where ( denotes the porosity and p denotes the interstitial fluid pressure. 
Note also that A and B denote different modes within the same solid. When the 
porosity of the elastic media is zero, the equation reduces to the fluid loaded elastic 
solid bi-orthogonality relationship. 

To illustrate how the bi-orthogonality property can be applied, we will examine 
the simple case of two fluids, each with different densities in contact at an interface 
along the y axis. In this case, the bi-orthogonality relationship reduces to standard 


orthogonality between sines and cosines. A diagram is shown in Figure(1), 





Figure 1. Bi-orthogonality of Two Fluid Layers 


where fluid 1 satisfies homogeneous Neumann|Ref. 8] conditions, p, = 0, and fluid 2 
satisfies homogeneous Dirichlet conditions, p = 0, along the lateral sides of the layers. 
From [Ref. 8] we know the eigenfunctions are sines for the Dirichlet conditions and 
cosines for the Neumann conditions. We equate the pressure and flow normal to the 


interface to obtain 














] 2 2 
2 S> (In + Rn)p\(y) | 7. | Timp’ (y) | -1? | qua) 
— Kid) - x?) 7 
U 7 | (In — Rn) UM (y) m | Tm2U2)(y) U 
where p)(y) and U)(y) are defined as follows: 
p(y) = cos(ny) and p(y) = sin(my) (11.28) 
UM(y) =cos(ny) and UL(y) = sin(my) 
Applying the bi-orthogonality condition as follows 
I U}) [p' = p*|dy 
[out =U" lay 
oe (11.26) 
[ U;" |p = p’|dy 


eS) 


pr[U = U"|dy, 


0 


we obtain the following two sets of equations, written in matrix form for easier ma- 


nipulation 


(11.27) 


We now define G = (Q7'G) and K = (J7~'K), which, if both formulations are to give 


identical answers, results in an identity GX = [. Solving for G and K we have: 


| * cos(py)sin(ny)dy 
Gon = (Q-'G) rate = 
| sin*(ny)dy 
0 
(11.28) 
| * cos(qy)sin(ny)dy 
| CP ee agg Oe a 
[ cos’(qy) dy 
leading to 
Gg = i a where p+ 7 must be even 
© ie Dy , 
(11.29) 
ay = ae eared where g+n must be even, 
which gives the matrix element (GX), as 
0 if p+q is odd 
> GpnKng = § 8(2 —5y1) n? 
n SS if 
7 Loni —@-ny NPT iseven 
(11.30) 


It can be shown that the above summation gives 


10 


(GH )pg = Opg (11.31) 


Knowing these relationships and given initial incident wave amplitudes, I, we can 
solve for the reflected and transmitted amplitudes from either one of the pairs of 


vector equations given by Equations (11.27). 


C. NUMERICAL METHODS 
iz Newton’s Method 


Newton’s method, one of the numerical methods used in this thesis, is a widely 
used method for solving nonlinear equations. We will briefly review the concepts and 
give a short example of this method. Newton’s Method uses a line tangent to the 
curve of a function and is based on the linear approximation of the function. We start 
with an initial guess zo, that we believe to be close to the root, and determine f (zo). 
Then we move along the tangent line until the point intersects the z-axis. This value 
of x now becomes our next iterate. In general terms, Newton’s Method is given as 


[Ref. 4] 





Oe) 
Inti =f F'(@n) ( ) 
where n = 0,1,2,.... Newton’s Method converges rapidly when the initial guess is in 


the neighborhood of the root. In fact, Newton’s Method is quadratically convergent. 


Here is an example [Ref. 4]: 


U (2) oop sinaeene: 
f(x) = 3+ cosx — e* 
We choose xp = 0 and perform 3 iterations as follows 
f (Zo) 
f (79) 


“0. = 36017 
i (zr) 


z 
oe — te nea = .3604217 


= .33333 





ZL, = Zo — 





Lg = ©) — 





J 


which is accurate to 7 significant digits. Although the convergence is rapid, there are 
other considerations that need to be addressed with this method. Newton’s Method 
requires 2 function evaluations per step. However, the total number of function 
evaluations is generally the same or less than other methods and Newton’s Method 
usually takes fewer iterations to determine a root with high accuracy. The main 
drawback to Newton’s Method is that it will not converge for certain functions. This 
happens when the tangent line is approximately horizontal or when we repeat a value 
previously used, for example x; = x4. Although these are strong considerations, 
this happens infrequently. For our needs, we will use this method to refine our wave 
speeds found in Chapter III. We will also use this method in refining our eigenvalue 


estimates found using the finite difference method which is discussed in Chapter IV. 


2. Finite Difference Method 

The finite difference method is used in numerically approximating derivatives 
and integrals of functions. We can use this ability to replace derivatives with difference 
quotients and thereby approximate solutions to difficult differential equations. There 
are three basic types of difference quotients; backward, central, and forward, with 
central being the most accurate. In certain instances, such as endpoints, backward 
or forward differences are often used due to the lack of additional points beyond 
the extremity of an interval necessary for accomplishing the central difference. In 
any case, it is customary to ensure that all difference operators have the same or like 
order of accuracy to avoid extra work and increase the possibility of error cancellation. 
Difference operators can also be applied more than once in order to obtain second 
order and higher difference approximations to higher order derivatives. For this thesis, 
we will concern ourselves with only second order differences. The central difference 


operators [Ref. 4] are as follows: 


12 


dx Ly41 — 27-1 


Fa | A 

(11.33) 
d* i41 — 22; + 24 
<7 i=t, = os de O(h?) 


where h = Az is the stepsize. Since forward and backward differences are O(h), we 
need to expand 2;4; and z;_; as follows to obtain formulas of higher order for our 


backward and forward difference operators: 


tig =U tah+aclh? + ac'ho4+--- 


(11.34) 
jy = 2; —-vh+c%h? — che +.-- 
Re-arranging terms we obtain 
2 — hy h h? 
ae ae ae ~zx, — x" + O(h) forward difference 
ae . - (11.35) 
zi = — ; — + ati - re + O(h?) backward difference 


WN 
2 


We will be able to substitute for the terms z7’, and x," in our final discretization. This 


will be discussed in more detail in Chapter IV. 
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IIf. DISPERSION EQUATIONS FOR 
ELASTIC HALF SPACES 


In this chapter we will examine several elastic fluid problems with infinite 
and finite layers. In each problem a relationship known as the Dispersion Equation 
is found. A system is dispersive if the phase velocity of a wave depends upon the 
frequency of the wave. Solving this dispersion relationship at a prescribed frequency, 
we can determine the wavespeeds for propagating modes and use this for an analysis 
of the system. These modes are also the necessary constituents in applying the bi- 


orthogonality relationships being studied. 


A. INFINITE ELASTIC HALF SPACE WITH A FREE 
SURFACE 


Free Space 


Infinite Elastic 
Half Space 
Figure 2. Infinite Free-Elastic Half Space 


The case of an infinite elastic half space with a free surface is described by 
Achenbach [Ref. 7]. We consider the two dimensional plane strain case. We apply 


stress free boundary conditions at the free surface: 


Toy = (Uy + vz) = 0 (11.1) 


Ty ee eee, — 0) iTe2) 


15 


Dividing Equation (III.2) by p and substituting our longitudinal and transverse wave 


speeds, c? and c} we obtain 


ee Myre ee Vii — Ol (irs) 


Assuming a solution for our vector potentials, y and y, as follows, 


= t(wt—kr) 
P= Pye” (111.4) 
y= w(y)e- i(wt—kz) 
our solutions for u and v become 
= 6, + We = thee 
re y a ia (11.5) 
UV = by — bz = Hy — thy 
and the two boundary conditions can be written: 
Uru = 0 
7 (111.6) 


(ci —2c4h)iku +c? v, =0 


Next, our potential functions must solve Helmholtz equations in the elastic solid, 


Veep tkiv =0 
V7 + ken =0 


(111.7) 


Dropping the tilde notation on the “amplitude” functions of y and the e~‘(“*~**) de- 


pendence, these functions satisfy the wave equation if 


2 0 
—— (II1.8) 
yp _ 7 aaa 0, 
where 7” = k* — k? and €? = k? — k4. This leads to exponential solutions 
ave eC foe 
a (111.9) 


y = Ce + De&Y 
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where we stipulate that Re(£,7) > 0. To ensure boundedness of the solutions as 


y —>+ —oo, we neglect solutions that “blow up” to obtain 


C= pe 
a, (III.10) 
== ie 


Substituting into the boundary conditions we have 


iknyp + 6? + tk(np — tk) = 0 
(cf, — 2h )ikliky + &B) + B(nPy — shed) =0 
Evaluating the amplitude functions at y = 0 and simplifying leads to 


(111.11) 


2k? — k2)\B = 20kED = 0 
( a (111.12) 
2iknB + (2k? —k#)D =0 
The determinant of the coefficient matrix in Equation (III.12) leads to the Dispersion 


Equation 


(Qh ks) atk En = 0 (111.13) 


In this instance, the waves are nondispersive because the phase speed determined 
by the above relationship is independent of the frequency. The frequency can be 
scaled out by noting k = . and dividing the equation by w*. Using cz = 5690m/s 
and cr = 3145m/s, the compressional and shear wave speeds for steel, we can apply 


Newton’s Method to solve for cr, the Rayleigh wave speed, and find 
CR = 2907m/s 


where cr < cr < cy. From [Ref. 7] we have the following approximate formula for 


computing the Rayleigh wave speed. 


s 862 + 1.14v 


~ III. 
CR pee (IIT.14) 
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Using this formula, we predict a Rayleigh wave speed of 2902.25m/s which is very 
close to our result. Now letting D = 1 we find that B = .6748:. At the frequency 
w = 1 and setting zc = y = 0, the real parts of the particle displacements (u, v) 
have been plotted in Figure(3). Note that this is the time history of the particle 
displacement at the surface of the elastic half space. In addition, as depth increases, 
the amplitude of the Rayleigh wave decreases and the motion alternates between 
retrograde and prograde. Rayleigh waves for various elastic solids appear in Table (1) 


later in this chapter. 


¥ Rayleigh Surface 
Vaves 


Figure 3. Rayleigh Wave Propagation 


B. FLUID LOADED INFINITE ELASTIC HALF SPACE 
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Infinite Elastic 
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Figure 4. Infinite Fluid-Elastic Half Space 


When a fluid is added, we have slightly different boundary conditions. Equa- 
tion (III.1) still holds, but Equation (III.2) is modified due to the fluid pressure at 
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the boundary. This leads to 


Tyy = A(Ue + vy) + 20, = —p (111.15) 


for which 


p= —p;s ot (111.16) 


and where ¢ is the acoustic velocity potential for the fluid and p; is the fluid density. 


The normal stress boundary condition becomes 


equ, + (ce = 2a, ue + iwi =) QUIN) 


We allow slippage along the fluid/solid interface (y = 0), but demand continuity of 


the vertical velocity there, i.e. Vsoia = Vfluia- This leads to a third boundary condition 


— iwy, — wk + dy =0 (111.18) 


The velocity potential function must also satisfy a Helmholtz equation 


pec a0 (111.19) 


Ww 
where ke = — and cp is the fluid wave speed. As before we determine the form of 
CF 


our solution using the wave equation and assuming a potential of the form 


d= dye et *) (111.20) 


Dropping again the tilde notation and the e~*“*-**) dependence, the amplitude func- 


tion must satisfy 


yy — Cd = 0 (111.21) 


where ¢* = k? — k2.. The amplitude function has the exponential solution 


Le 


O= Peumimue (111.22) 


Disregarding solutions that “blow up” as y — too we obtain the amplitude functions 


for the solid and fluid 


O— peu 
b= Deb¥ (Ihe 35 
d= Ee7%¥ 


Substituting these into the three boundary conditions we have at y = 0: 


2ukny + (€? + k*)yb = 0 
(ch? — k(c3, — 2ch))yp — 2ikEcReb + iw8Lg = 0 (111.24) 
ny +wkp —Coe=0 
Simplifying we obtain 


2ikn B + (2k? — k2)D = 0 
(2k*c} — ki ch) B — ukéceD + wt E = 0 (111.25) 
w7B+wkD—-CE=0 
Computing the determinant of the matrix and simplifying results in the Dispersion 
Equation. However, unlike the previous case, we now have an added term due to the 


presence of the fluid half space. This equation is 


(2k? — k2.)? — 4k? né] + rare = (111.26) 


Ss 


Note that once again the frequency can be factored out of the above expression which 
implies a nondispersive media. Applying Newton’s Method using the wave speeds for 


steel we determine what is referred to as the Scholte[Ref. 7] wave speed, 


c, = 2903.567 — 35.3252 m/s 
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Note that in the case with the fluid half space we have a complex wave speed. Letting 
E = 1 we find B = 2.51 + .3952 and D = —.030 + 3.892. Setting w = 1 and plotting 
the real part of the horizontal and vertical displacements at z = y = 0, we obtain the 
particle displacement history as this Scholte wave passes through the point z = y = 0. 


The plot looks the same as Figure (3). 


C. FINITE FLUID INFINITE ELASTIC HALF SPACE 


Fluid 


Infinite Elastic 
Half Space 


Figure 5. Fimite Fluid Infinite Elastic Half Space 


We next consider a fluid layer of finite thickness overlying an elastic half space. 
The boundary conditions from before, Equations (III.15), (III.17), and (III.18) still 
hold. However, an additional boundary condition is needed at the boundary between 
the fluid and free space which is that the pressure be zero there. From Equation 


(III.16), this leads to 
p= aoa — lg) CieZ 7) 
Also, from before we obtain the solutions to the Helmholtz equation where 


yy = Ae” 
y= Bef 


(111.28) 


Our third potential changes because we now have a finite layer. We are unable to 


discard one of the solutions, but instead must have 


2] 


= Ce + De“ (111.29) 


Simplifying this by combining the terms, we choose 


$ = Csinh(¢(h — y)) (111.30) 


Notice that by making this choice we have incorporated our additional boundary 
condition so that when y = h, 6 = 0. Substituting into the boundary conditions at 


y = 0, we have 


2ikny + (€? + k*?)b = 0 
(cin? — k* (ct — 2ch))y — 2tkécep + wo = 0 Gieet 
wryyp + wkp — CC cosh(C(h — y)) = 0 


Or 


2iknA + (2k? + k2)B =0 
(2k?c%. — w?)A — ike B+ tw sinh((h)C = 0 (III.32) 
wnA +wkB—Ccosh(Ch)C = 0 
Computing the determinant and simplifying leads us again to a Dispersion Equation, 


which is 


(2k? — BB)? — 4k ne] + ran tanh(Ch) = 0 (111.33) 


Unlike the previous case, the frequency can no longer be factored out. This dispersion 
relation signals that we now have a dispersive media which is intermediate between 
the other two limiting cases. Notice that as w — 0, the “fluid term” is O(w°) while the 
solid terms are O(w*), and therefore the wave acts like a Rayleigh wave. As w — oo, 
tanh(¢h) — 1 we obtain the Scholte limit. This makes sense physically because as 
w—+0,A~ O(=) becomes large and the fluid layer will appear to become extremely 


“thin” relative to the propagating wavelength. On the other hand if w — oo, the 


Zz 


fluid will appear to be “thicker” in terms of the number of wavelengths that can fit 
onto the fluid layer, and we get the other limit of an infinite fluid layer. As in the 
previous problems we solve for the Scholte interface wave speed. A height of 1m is 


used for the fluid layer in the computation. The wave speed for steel is 
c= 2906 — .0011z2 m/s 


Letting C = 1 we find that A = .00947 — .104r107®2 and B = .147r1078 — .005711:. 
We now setw = 1,h=1, andz=y=0. This plot is the same as our Rayleigh wave 
plot in Figure (3), which was expected. Listed below in Table (1) are Rayleigh and 


Scholte wave speeds for various elastic solids. 


7 Steel 5690 | 3145 | 2007 | 2904 — 35; | 2006 — .001r 
[ron (Cast) [4175 [9308 [2133 | 2127 — 451 | 2133 — 00007 
[Aluminum | 6242 | 3144 | 2030 | 2904 — 1121 [ 2930 — .0031_ 
[Glass (Pyrex) || 5637 [ 3207 | 3026 | 2003 — 1530 | 3026 — 0047 
[Rubber (Hard) [27 | 864 | sid oor | sid 


Table I. Rayleigh and Scholte Wave Speeds for Various Elastic Solids 








D. FLUID LOADED ELASTIC PLATE 





Figure 6. Fluid Loaded Elastic Plate 


We consider an elastic layer of finite thickness. We examine the relationhip 


of a finite fluid media in contact with a finite elastic medium or plate. As shown 
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in Figure (6), the origin is assumed to go through the middle of the elastic solid. 
This will help later to show the symmetric and anti-symmetric modes as discussed by 
Achenbach [Ref. 7] in the unloaded, “free plate” problem. Our boundary conditions 
are slightly different with this case. We choose a thickness h, for the elastic layer and 
a thickness hz for the fluid layer. At y = = we have 


Try = 0 


(111.34) 


Bye L 


h 
AL y= ae we have 


av 


Tyy = —p (IIT.35) 


Usolid = Uf luid 


Finally, at y = ho + 2 we want the pressure to again be equal to zero. We again 
assume solutions to the wave equation. However, we now have a finite fluid and elastic 
layer, and therefore must retain all of the exponential terms. We write the solutions 
in terms of cosh and sinh in order to better see the symmetric and anti-symmetric 


portions of the solutions. We have, 


(ny) + B sinh(ny) 
w = Ccosh(€y) + Dsinh(€y) (111.36) 


¢ = Esinh(¢(h — y)) 


yy = Acosh 


h ; ae 
where h = — + hg, and we have taken into account the boundary condition that 
D—Wrateie- 


Substituting into our boundary conditions and letting p = tii we have 
Ps 
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2iknyy(—4) + (2k? — kp) y(-4) = 
2iknyy() + (2k? — kp)d(F) = 
(2k? — kR)yp(—) — kev, (-#) = 


I 


(2k? — kp) o(B) — BWhby (GB) + AP b(hs 
iwnyy(Z Ver wkp( 2) — dy(h2 
Trying to find and simplify the determinant is difficult. We try a simpler approach by 


a) 
0) 
ce UU 
eal 
a) 


adding and subtracting like equations to obtain our Dispersion Equation. Adding and 
subtracting the first two equations and doing the same with the next two equations, 


we have the following 


2:knB cosh(n*) + (2k? — k%)C cosh(€2) = 
2iknA sinh(n2) + (2k? — kp) D sinh(€2) 
(2k? — kp) A cosh(n4t) — 2ik€D cosh(€4) + $22 E sinh(Che) 
(¢h2) 


hy 
2 

hy 
2 


| 
a eee 


ba (111.38) 
(2k? — ky) B sinh(n4t) — 2ikEC sinh(€41) + $22 E sinh(Che 
iw yy(F) + wkb(F) — gy(ho 


From the first two equations, we solve for A and B in terms of C and D. Substituting 
into the next two equations we obtain a relationship of C and D and find F in terms 
of C. This is not as difficult as it seems, but does require accurate bookkeeping to 
keep up with the algebra. Substituting all of these into the last equation, and after 


much algebra, we obtain the Dispersion Equation for a fluid loaded plate. 


[(2k? — kp)? — 4k*né tanh(n*) coth(é4)][(2k? — k7)? — 4k?n€ coth(n*) tanh(é*)] 
= =eckbn anh( Cha) a — ki)? coth(nh,) — 4k?n€ coth(€h )] 
& (111.39) 
Because the dispersion relationships will involve infinitely many complex roots, due 


to the finite nature of the elastic/fluid system, Newton’s method may indeed be able 


to find one or several of the roots, but will not be efficient in finding the first N roots 
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(ordered by magnitude), simply because required starting values are not available. 
One lacks the confidence that all of the eigenvalues will be found with Newton’s 


method without sufficiently close starting values. 
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vn RESEARCH RESULTS 


A. DISCRETIZATION PROCESS 

We begin the numerical analysis of the fluid loaded plate by applying finite 
difference approximations to the boundary value problem. We know that in the fluid 
and solid, the respective reduced Helmholtz equations must be satisfied. We define 
(N + 1)hy as the thickness of the fluid layer and (M + 1)h, as the thickness for the 
elastic layer, where (V + 1) is the number of discrete points in the fluid layer and 
(M + 1) as the number of grid points in the elastic layer, and hy and h, are the 
stepsize for the fluid and elastic layers, respectively. In order to obtain our unknown 
modal wave number, k?, on the right hand side of the system of equations, several 
substitutions are necessary. We define /* = wz, and also substitute ¢; = ¢. Therefore 
in the final analysis we will be solving for normal modes of the functions ¢;,y, and 


w*. Dropping the superscripts and subscripts to avoid clutter in the discretization, 


we have 


bf + ktd: = kid; w= 1..N —] 
ol + kiy; = k2y; = Vea — | (Tiah) 
yl + kab; = kev; 2=1..M—1 
The additional boundary conditions on the fluid loaded elastic plate from Chapter 
III, can be reduced after completing the above substitutions and applying the finite 


difference method, to the following: 


on = 0 

20" + (2kt — kp )yp — 2p! = Lo 
—2k2y! + kaa + 2p” =0 
b-y=s 

Quo" + (2k? — k2)p — 2a’ = 0 
—2k2p' + ko + 2” = 0 


(IV.2) 
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where z = 0 is the interface between the fluid and the solid. In order to incorporate 
these additional boundary conditions into Equations (IV.1), we must solve for what 
we will call, pseudo nodes. A pseudo node is a point outside of the fluid or elastic 
solid that will be used in the differencing of a potential. This can be seen below in 


Figure (7), 





Figure 7. Discretization of Fluid Loaded Elastic Plate 


where ¢; is the pseudo node for the fluid and y_1,%_1, ¢yg41, and Wy are the 
pseudo nodes for the elastic solid. 
The equations for the pseudo nodes require some careful algebra and can be 


reduced to the following: 





oe as ae: v we — 2hyw*dy + = [1 + he (kp — k7)]y0 — : shsk?w*y, — te po 
= [1+ A3(kE — kz) ]~0 + (1 — AZkz) pi — G2 + (oh, + EPS obo — hey + Bet ec 

Waa = (44 kN?) vo — 3di + (Aske — FZ) yo + (GE — WkEhs)o1 — Gye + bps 

pm+i = [1 — h2(k? — ke) om + (1 — h2k7)em-1 — em-2 + Dhavaer it (Gas — 2hs)vm 


vm+i = (4— Akt) bm — 80m-1 4+ (~ — Ask) ym + (2k7Zhs — =) eM—-1 + EH PM-2 
(IV.3) 


These pseudo nodes are substituted into our Helmholtz equations from (IV.1) to 
obtain equations for the boundaries 2 = 0 and 2 = M. Our discretized boundary 
equations are now entered into a MATLAB code given in the Appendix. Note that 
for the MATLAB code we must start our index at 1 as MATLAB does not have 0 
indices. A shift is required so that the first N places correspond to ¢, the next M+1 
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places correspond to y, and lastly, the next 1f+1 places correspond to =. This allows 
us to keep track of each of the three potentials. MATLAB computes our eigenvalues 
and orders them in ascending magnitude. The code also calculates the normal and 
shear stresses for use in applying the bi-orthogonality condition. Finally, various 
graphs involving our discretized solutions will be plotted against analytical solutions. 
The analytical solutions for stresses, pressures, and displacements are determined 


using the coefficients listed in the last section of Chapter III. 


B. NUMERICAL RESULTS 
1 Modal Wave Numbers 


The elastic solid used for the following analysis will be steel. We begin with 
an examination of the eigenvalues of our discretized matrix. In the discretization, 
w = 2r, N= 10, and M = 50, leading to a square matrix of size 112 from which 112 
eigenvalues are determined. Since we are solving for k?, we must take the square root 
of the eigenvalues and place them in order. Determination of the accuracy of these 
wave numbers will involve substitution into Equation (III.39), our fluid loaded elastic 
plate dispersion relation. A plot of the norm of the eigenvalues is shown in Figure 
(8) on the following page. Aside from the last few eigenvalues, we get a nice ordering 
of the eigenvalues. 

We attempt to refine our approximate wave numbers using several iterations of 
MAPLE coded Newton’s Method with the eigenvalues as initial iterates. The results 
are shown in Table (II). For some of our “refined” wave speeds, Newton’s Method 
attempts to converge to a root. Newton’s Method diverged for the roots marked with 
an asterik. After refinement, several of the roots show improvement in satisfying this 
equation, 1.e. the magnitude of the residual of the relationship is much closer to zero. 
The fundamental root has a residual magnitude change from 10~7° to 107°°. The 
second and third roots change almost as much. However, the rest of the roots exhibit 


less of a change in magnitude, indicating that the starting values for the roots may 
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Norm of Eigenvalues 





0 0.05 0.1 Os 0.2 0.25 0.3 0235 0.4 
magnitude of lambda 


Figure 8. Eigenvalues of Discretized Matrix 
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not be close enough to the actual roots for Newton’s Method to converge. 


Refined k, 


0.2804 
Pa [oer — 08003 | 
75 [006711 + 008003" |_| 
8 [ .004703 — 01666" | iY 
C7 [004703 +.01666" | 
Ts 002550 — 031i | ola 02i0 | 4855 
9 [0025504 031i | 01127 +0210; | 4355 —_+| 


Table II. Calculated versus Refined Modal Wave Speeds 


















Table (II) shows the first three roots converging to a real root.The other “refined” 
k,,’s do not appear to converge to any of the k,,’s used as initial starting values. It 
may be necessary to search more than the first ten modes to help determine if a shift 
of the “refined” values to a corresponding k, is possible. For this reason, we will 


concentrate on the first few modes. 


2. Stress Potentials 

The following graphs show the analytical versus the discretized solutions for 
the fluid and solid potentials ¢, y, and w. Figure (9), displays the best results. The 
analytical solution of ¢ is almost identical with the discretized solution. Note that 
the discretization code has ¢, as the first interior point and not the boundary. This 
accounts for the slight mismatch at the fluid-free space boundary. Figure (10) shows 
a simliar pattern where our solution is nearly identical, but diverges slightly as we 
near the fluid-elastic boundary. Figure (11) diverges at both ends but is on a scale of 
10-°. When viewed on a small scale, the discretized solution does diverge more than 
gm and ». It is obvious that there is a problem at the fluid-elastic boundary and from 


Figure (11), there is a slight indication of a problem at the bottom boundary. 


3] 


analytical (0) —- discretizes (x) 


Analytical vs Discretized for Phi (fluid) 





Figure 9. Analytical versus Discretized Solutions for @ 
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Analytical vs Discretized for Varphi (solid) 
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Figure 10. Analytical versus Discretized Solutions for y 
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Analytical vs Discretized for Psi (solid) 
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Figure 11. Analytical versus Discretized Solutions for w 
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Figure (12) illustrates a pressure-stress graph. Notice that at the boundary 
of the fluid and elastic layer a discontinuity exists. One of our boundary conditions 


from the fluid loaded plate is 


Tyy = —P (IV.4) 


This means that our plot should connect smoothly at the interface. The discontinuity 
indicates a problem with the discretization at the interface. This problem was evident 
in the previous plots and is re-enforced with this graph. 

Figure (13) plots the wave numbers of the first three modes as a function of 
increasing frequency. We see that our fundamental mode is purely real and our second 


and third modes are complex. 


3.  Bi-orthogonality 
The bi-orthogonality condition of differing modes from Chapter II indicates 
that we should expect a diagonal matrix in forming the bi-orthogonality product of a 


set of modes, where the diagonal element, x, is the “norm” of each individual mode. 


0 hy 
In = a [ri yl) - Tey) tactic dy -| (pM U™ | piady (IV.5) 
—f2 


The bi-orthogonality condition states that the off-diagonal elements, the product of 
two separate modes, should go to zero. Numerically, we expect these elements to tend 
toward zero. The numerical bi-orthogonality relationship for the fluid loaded plate 


yields the following matrix: 


AAT OO ers hill a De eile hOO cal rmeneggo x 10° 
es lean eed i aE Oe xaos e 30x 10> 232639" < 110° 
ieee eel N02 46a 10° 2670 x 10°. 2.670 x 10° (IV.6) 
2 Olea OOR lO iG oO 10° ns24 x 10? 1.832 x 10° 
Uj Ueno 0 = Gelso 10°) 1.832% 10° 1.524 x 102 
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Wave Propagation for Stee! 
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Figure 13. Wave Mode Propagation - Modes 1-3 


on 


Notice in this matrix, most of the off-diagonal elements tend to exhibit some form 
of reduction. However, this did not happen as quickly as believed. This indicates 
that the bi-orthogonality condition is very sensitive to small errors in calculating the 
modes. A more careful numerical treatment might improve the diagonal form of the 
matrix, but it may be necessary to limit the application of the numerical method 
to eigenmode determination only, from which the eigenmodes can be refined and 


ultimately found by analytical means. 
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Ve CONCLUSIONS AND FUTURE 
RESEARCH 


We can draw several conclusions on the discretization of the fluid loaded elastic 
plate. The eigenvalues, or modal wave numbers, that resulted from the discretization 
nearly satisfy the Dispersion Equation (III.39). Using Newton’s Method, we deter- 
mined that our initial values may not be close enough in some instances to the actual 
roots as is evidenced from the lack of improvement of the residual calculation in the 
“refined” root. Therefore, we conclude that there are numerous other roots that we 
have not found. 

The biggest disparity lies in the plot of pressure and normal stress versus the 
height of the fluid and elastic layer. As seen in the discretization code in the Ap- 
pendix, differencing of the potentials is necessary when solving for the stresses and 
displacements. Therefore, any numerical error in the fluid-elastic boundary differenc- 
ing is compounded. Recall, in theory the normal stress is equal to minus the pressure 
at this boundary. The discretization does not show this. Finally, we conclude that 
our graph of the first three modes in Figure (13), could indeed be modes for this 
problem. However, they are most likely not the first three propagating modes. 

The trouble with the results is mostly due to the following: a) there are errors 
in the discretization or b) the discretization is unable to properly model the boundary 
conditions. Indications of the first problem are evidenced in two ways. The first is the 
inability of the eigenvalues to converge to the correct analytical eigenvalues. Secondly, 
the plotting of the corresponding eigenvectors do not satisfy the boundary conditions 
as required. The solution to this problem is to find the correct way to difference the 
boundary conditions. Nothing can be done to fix the second problem, except to try 
a different discretization. The code as it stands does find some roots, but until the 
boundary conditions can be satisfied, there is no real confidence that the roots are 


comprehensive. 
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Several new problems arise from this research. The first is to determine a 
better way of differencing the boundary of the fluid-elastic layer. Once this problem 
is corrected, we can solve for the coefhcients of transmitted and reflected waves given 
and incident wave. This is useful when extending the elastic case to the more complex 
porous case. If the porous case can be solved numerically, this will allow us to create 
ocean acoustic models where the transmitted and reflected waves can be determined 
from known coefficients. This has potential applications in sonar and mapping of 
the ocean floor. Finally, a similiar problem involves modeling a fluid layer only, with 
different sound velocity profiles and an acoustically hard boundary. The discretization 
would not involve the elastic boundary and may be a better way to approach the fluid- 


elastic case. 
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APPENDIX. MISCELLANEOUS DATA 


1 MATLAB CODE FOR GENERATING WAVE SPEEDS 


functions Ll) =~elastic(E, nu, rho) 


YC GC ESE AI IIA a aE a a ak 
7, ELASTIC - This function computes the elastic constants and 


i wavespeeds for various solids. Input values are the 
hh Young’s Modulus, Poisson’s ratio, and the density. 
yf 


7, Copyright 1997 by C. R. Myers, III 
YSERA ICRA IEE 


lambda=0O; 
mu=0O ; 
cl=0; 
ct=0; 
E=E*10710; 


lambda = (E*nu)/((itnu)*(1-2*nu) ) 


mu = E/(2*(i+nu) ) 
cl = round(sqrt((lambda + 2*mu)/rho)); 
ct = round(sqrt(mu/rho)) ; 


2. DISCRETIZATION MATLAB CODE 


(METI TTPO CTOCSCCSSOCS CSCC CLOPSCCOCL OCLC SOCS OCLC STP SSS LT ELS eee: 
4, Fluid-Elasic Finite Difference Matlab Code 


4 Copyright 1998 by Coley R. Myers, III 
Y, GAGA GARG GRR RA IER aK a a a ok ak 2 ok 


Aeinitialize 


clear; 

format long e; 
w=0; 

E_val=[]; 


4 allows program to run in loop for obtaining 
4 modes over a specified frequency range. 


4] 


counter=1; 


4 loop is a counter and is used to obtain graph 
4 graph of modes over a specified frequency range 


for loop=1:counter 


/, Variables input for specific elastic solids 
Y, COO KG aR a a Rak ik kkk a ak ak kk 


Seenq nial) 3 

Lf=10; 

Ls=200; 

w=2*pit(w+.1); rf=1026; rs=7700; 
ct=3145; cl=5690; cf=1500; 
modes=10; 

mu=rs*ct” 2; 

lam=rs*cl°2 - 2*mu; 

hf=1; 

hs=4; 

N=Lf/hf ; 

M=Ls/hs; 

d=N+M+M+2; 

A=zeros(d); 

kf=w/cf; kt=w/ct; kl=w/cl; 


4% Ensures Matrices are clear when loops 1 


Vi=() ;Phi=(] ;Vphi=(J] ;Psi=(] ;Psixy=() ;Vphiy=(J ; 
v=] ;D=0] xPres=[] ;u=(); 

u=(] ; v=(] ;Styy=(] ;Stxx=0] ;Stxy=(01; 

pr_styy=(] ;f1=0 ;F1=0) ;Y=0;y1=0 ; 

Vphia=[] ;Psia=D) ;Phia=(] ; 


7% Fluid BC Phi(0)=0 and 1st equation 


p=1; 
NCon i) ==(2/ht 2-kio 2): 
A(p,2)=1/hf°2; 


p=pti; 


for 1=2:N-1 
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ACp,i-1)=1/hf~2; 
Alp: i) == (2/ht«2-kr 2) ; 
Atp,it1)=1/h£t 2; 


p=pti; 
end; 
4, 3rd BC for Phi 


A(p,N-1)=2/hf7~2; 

A(p ,N)=-2/hf*2+kf*2+ (hs*w* 2*rf ) /(mu*hf ) ; 
A(p,N+1)=(w72/ (hs*hf))*(1+hs*2*(kt*2-k1°2)); 
A(p ,N+2)=-hs#w* 2*k1* 2/hf ; 

A(p,N+3)=-w72/ Chs*hf) ; 

A(p ,N+M+2)=(kt~2*hs* 2*w72)/(2*hf) ; 

ACp ,N+M+3)=-2*w7 2/hf ; 


4 1st BC for Varphi 


p=pti; 

A(p ,N+M+2)=(2/hs)+kt“2*hs/2; 
A(p ,N+M+3)=-2/hs; 
A(p,N+3)=-1/hs*2; 

A(p ,N+2) =(2/hs*2)-k1°2; 
A(p,N+1)=kt*2-1/hs*2; 

A(p ,N)=rf/mu; 

p=pti; 


% 2nd BC for Varphi 

for i=(N+2): (N+M) 
A(p,i-1)=1/hs72; 
Gog 2 oe 2a 2) 


A(p,it+1)=1/hs*2; 
p=pti; 


end; 


“4 3rd BC for Varphi 
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A(p ,N+2*M+2)=(kt*2*hs/2)-2/hs; 
A(p,N+2+*M+1)=2/hs; 

A(p ,N+M+1)=kt~2-1/hs*2; 
A(p,N+M)=2/hs*2-k1°2; 
A(p,N+M-1)=-1/hs“2; 

p=pti; 


% ist BC for Psi 


A(p,N)=rf/(nuths) ; 

A(p ,N+1)=kt*2/hs-2/hs“3; 
A(p,N+2)=4/hs*3-2*k1*2/hs; 
A(p,N+3)=-2/hs°3; 
A(p,N+M+2)=2/hs*2; 

A(p ,N+M+3)=-2/hs*2; 

p=pti; 


7% 2nd BC for Psi 
for i= (N+M+3) : (N+2*M+1) 


A(p,i-1)=1/hs*2; 
A(p,i)=-(2/hs*2-kt~2) ; 
A(p,i+1)=1/hs*2; 
p=p+1; 


end; 
VA betel ()2(; aeons 12st 


A(p,N+M-1)=2/hs°3; 

A(p, N4+M)=(2*k1*2/hs)-4/hs“*3; 
A(p,N+M+1)=(2/hs*3)-kt*2/hs; 
A(p,N+2*M+1)=-2/hs*2; 
A(p,N+2*M+2)=2/hs*2; 


7 Final matrix is calculated 
[V, D]=eig(A); 


D=sqrt (diag(D)); 
eet! |=sont op). 
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7%, Lambda holds the number of modes specified for 
% Uusemwith later calewulations 


Lambda=Y(1i:modes) ; 


4 when loop > 1, E_val holds values for graph 
/, of modes over a specified frequency range 


E_val = [E_val Lambda] ; 
7, Orders eigenvectors 


for 1 = 1:modes 
eee’: 
WES al WC gap ial s 
end 


4, Separates stress potentials for calculations 


for 1 = 1:modes 
Php ea eV eC taNere a 
Vphi = [Vphi V1(N+1:N+M+1,1)]; 
Psi = [Psi V1(N+M+2:d,i)]; 
end; 


yA 2 2 2K 2K 2K 2K 2K 2K 2K 2K 2K 2K OK OK 2K OK OK 2K 2K 2K 2K 2K 2K OK 2 2K 2K OK OK 2K OK 2K 2K OK RK OK 2K OK 2 2 2 2 2K 2K OK OK OK 2K OK OK OK OK OE OK OK OK OE OK OK OK OK 


7, Solve for Pressure 


Pres = -rf.*Phi; 
U = -(1/(CI*w)) .*Phi; 


if 2K 2K 2K 2K 2K OK OK 2 2K OK 2K OK 2K OK KK 2K 2K OK 2K OK OK 2K 2K 2K OK 2K OK 2K 2K OK 2K OK OK OK OK OK 2K OK OK OK OK OK OK OK OK OK OK OK OK OK OK OK OK OK OK OOK KOK OK KK 


4 Differencing for discretized potentials 
eta2 = Lambda.”2 - k1°2; 

mi2 = Lambda. 2 — kt 2; 

zeta2 = Lambda.”2 - kf72; 

4 Finding Psi_xy (same as Psi_yx) 

j=l; 


for i=1:modes 
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Psixy(j,1)==-C(PsiQGjtiga) =) 4 ahs 272) +x12 Gee ae 
(hs + (hs°3/6) .¥xi2(1)))- 
end 


for j=2:M 
for 1=1:modes 
Psixy(j,i)=(Psi(j-1,i1)-Psi(j+1,i))/(2*hs) ; 
end 
end 


j=M+1; 

for i=1:modes 
Psixy(].1)=(Psi(j-1,1)- (1 + (hs° 2/2). *xi2(1)) -+Psi Gee 
(hs + (hs73/6) .*xi2(i)); 

end 


7, Finding Vphi_y 


H- 15 

for 1=1:modes 
Vphiy (j ,1)=-C(Vphi(j+1,i)- (1 + (hs*2/2) .*eta2(i)) .* 
Vphi(j,i))./Chs + (hs*3/6) .*eta2(i))); 

end 


ton j-2:M 
for 1=1:modes 
Vphiy(j,i)=(Vphi(j-1, i)-Vphi(j+1,i))/(2*hs) ; 
end 
end 


j=M+1; 

for 1=1:modes 
Vphiy (j,i)=(Vphi(j-1,i)- (1 + (hs72/2).*eta2(i)).* 
Vphiey i) )o/ (hs + (hs 3/6) .#eta2G)): 

end 


TOO OGG GRR IIR I Ka ak 2k ak 2k 2k ak ak 


4 Solving for Tau_yy, Tau_xy, and Tau_yy 
for 1=1:modes 


Styy = ([Styy (-((lam + 2*mu)*k1°2).*Vphi(:,i) - (2*mu).* 
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Psa yea): 


Stxx = [Stxx ((2*mu).*Psixy(:,i) - ((2*mu)*Lambda(i)“2 + 
lam¥k1~2) .*Vphi(:,i))]; 
Stxy = [Stxy (((2*mu*I)*Lambda(i).*Vphiy(:,1i)) + 
(2*Lambda(i)7~2 - kt72)*(1/(I*Lambda(i))).*Psi(€: ,i))]; 
end ; 


4 Solving for u and v 


for 1=1:modes 
u = [u (I.*Lambda(i) .*Vphi(: ,i)+(1./(1.*Lambda(i))).* 
Remax Gs sa) )t: 
wee iv (Vpniy(:,1) — Psi(;)4) ae 

end 


(U.K RR a a a EI a a a6 2k a ok a6 ok 
/, Applying bi-orthogonality condition 


for 1 = 1i:modes 
Pile="(Pres(:51)- Sty. 1) 1; 
pr_styy = [pristyy F1]; 

end 


£1 = (Stxx’*u - Stxy’*v) - Pres’+U; 


for j=1:modes 
for 1=1:modes 
Poca t Cen econy ch ici.) ))), 
end 
end 


i AE 2K 2 2K 2 2K 2K ok ok 2k 2k 2 ok 2k 2 2 2K 2K 2 ok KK 2k KK 2k 2K 2K 2K 2K ok 2 2 2 2 2 2 2K 2 2 2 2 2 2 2 2K 2 2 OE 2 2 OK 2 2K ok 2K oe ok 2K 2K ok 


/Z Plot of norm of eigenvalues 


mi=sqrt(Y.#conj(Y)) ; 

m= 562 = 121: 

figure 

plot (Y¥1(1:56),i,’o’) 
title(’Norm of Eigenvalues ’) 
xlabel(’magnitude of lambda’) 
ylabel (’modes’ ) 
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fA I A a a a a ak ak ak ak ak ak 2k ak ak 2k a ak 2k ak ak ak ake ak ake ake ak ak ak akc 2 ak 2k akc ak ok ak ak ak ak ak ak ak ok ok 2k a 2 2 


/, Determine coefficents for analytical solutions 


k=Lambda(1,1); 

Getaze=— ik 24—- kie2: 

ax 2 ve ke 2 = kt 22: 

dazetaze= k 2 kz 2; 

hs=200; 

hf=10; 

h=hfths; 

7, Change value for C 

C=.00004; 

D1i=(2*k*2-kt*2)*2*tanh (sqrt (deta2)*(hs/2)) ; 

D2=4*k* 2*sqrt (deta2) *sqrt (dxi2) *tanh(sqrt (dxi2)*(hs/2)) ; 
D3=(2*k*2-kt~2)“2*coth(sqrt (deta2) *(hs/2) )*tanh (sqrt (dxi2)* 
(hs/2)) ; 

D=(D1-D2) /(D3-4*k*2*sqrt (deta2) *sqrt (dxi2)); 
B1i=-(2*k*2-kt*2)/(2*1*k*sqrt (deta2)) ; 

B2=(cosh(sqrt (dxi2) *(hs/2))/cosh(sqrt (deta2)*(hs/2)))*C; 
B=B1*B2; 

A2=(sinh (sqrt (dxi2) *(hs/2))/sinh(sqrt (deta2)*(hs/2)))*D; 
A=B1i*A2; 

E1=(D1-D2)#C; 

E2=ct/(k*kt*sqrt (deta2)*(rf/rs) ) ; 

E3=cosh (sqrt (dxi2)*hs/2)/sinh(sqrt (dzeta2)*(hf+ths/2)); 
E=E1*E2*E3; 


Vesa = 217: 
9 / LES a ale Mie 
WON Salis 


Vphia = A.*cosh(sqrt (deta2.*y1)) + B.*sinh(sqrt (deta2.*y1)) ; 
Psia = C.¥*cosh(sqrt (dxi2.*y1)) + D.*sinh(sqrt(dxi2.*y1)) ; 
Phia = E.*sinh(sqrt (dzeta2) .*(N-y2)) ; 


YSIS GRR RE 
4, Plot analytical vs discretized solutions 


figure 


Vphia=sqrt (Vphia. *conj (Vphia)) ; 
ples Wphiasyi 07 6N pat hae 
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title(’Analytical vs Discretized for Varphi (solid) ’') 
eabeld2 x.) 
ylabel(’analytical (0) - discretizes (x)’) 


figure 

Psia = Psia.*I.*k; 

Psia = sqrt (Psia.*conj(Psia)) ; 
PlobG@Psia- vile Ommpsi@el evi,’ x”) 
title(’Analytical vs Discretized for Psi (solid)’) 
xlabel(’x’) 

ylabel(’analytical (0) - discretizes (x)’) 


figure 

Phia = I.*w.*Phia; 

Phia = sqrt (Phia.*conj(Phia) ) ; 
piltom@rhiae2. O° Phi (1) ,y2,'x") 
title(’Analytical vs Discretized for Phi (fluid)’) 
xlabel(’x’) 

ylabel(’analytical (0) - discretizes (x)’) 


yf 2K FE 2K OK 2 2K KE 2 IE 2K 2 IK OK OK eK 2K 2 eK 2 OK 2 2K 2g KE KE eK KK EK cE 2 oe KK 2g 2 cK 2 IE 2 EE 2K 2 KE 2K KK 2K ok 2k 


% end loop 


end 
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